function [ssr,out] = ObjFun_cMSM(params,knotvec_c,student_id_l,Q,T,ThetaE,contract,Q80freq,pay,SimTimeCum,Activeids,Marginalids1,Marginalids2,...
                                Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,betaq1,betaq2,betaq3,knotvec_q,Display,bootstrap) 


%% 

%%%%NOTES:
%%%%bootstraps(1)<10 means pattern search algorithm in use (as opposed to genetic algorithm)
%%%%bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11 means point estimates (as opposed to bootstraps)
%%
K_c          = length(knotvec_c) - 1 ; %#ok
temp         = size(params) ;
if temp(2) == 1
    pi_c = params ;
else   
    pi_c = params' ;
end

%%%%WE USE THE CONDITIONAL DISTRIBUTIONS OF POSITIVE WORK (I.E.,
%%%%LEARNING TASK ACCOMPLISHMENT >=1) FOR THE CONTRACT-3-EQUIVALENT SETS OF STUDENTS
%%%%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%% IMPORTANT NOTE: THE SIMTIMECUM MATRIX IS ARRANGED SO THAT THE FIRST LENGTH(ACTIVEIDS) SHEETS
%%%% ARE FOR PEOPLE WHO COMPLETED 2 TASKS, THE NEXT LENGTH(MARGINALIDS1) SHEETS ARE FOR STUDENTS IN
%%%% CONTRACT 1 WHO COMPLETED EXACTLY 1 TASK, AND THE LAST LENGTH(MARGINALIDS2) SHEETS ARE FOR
%%%% STUDENTS IN CONTRACT 2 WHO COMPLETED EXACTLY 1 TASK.  THUS, WE START HERE BY ESTABLISHING THE
%%%% SAME ORDERING ON Q, T, BASE, WAGE, AND CONTRACT (AND DISCARDING ALL OTHER ROWS OF THOSE 5 ARRAYS):
%%%%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
NActive          = length(Activeids) ; 
NMarginal1       = length(Marginalids1) ; 
NMarginal2       = length(Marginalids2) ; 
student_id_ltemp = nan(NActive+NMarginal1+NMarginal2,1) ; % = student_id_l( (Q>=2&contract==3) | (Q>=1&contract==2) |(Q>=1&contract==1) ) ; 
Qtemp            = nan(NActive+NMarginal1+NMarginal2,1) ; %  Q( (Q>=2&contract==3) | (Q>=1&contract==2) | (Q>=1&contract==1) ) ; 
Ttemp            = nan(NActive+NMarginal1+NMarginal2,1) ; %  T( (Q>=2&contract==3) | (Q>=1&contract==2) | (Q>=1&contract==1) ) ; 
ThetaEtemp       = nan(NActive+NMarginal1+NMarginal2,1) ; %  ThetaE( (Q>=2&contract==3) | (Q>=1&contract==2) | (Q>=1&contract==1) ) ; 
contracttemp     = nan(NActive+NMarginal1+NMarginal2,1) ; %  contract( (Q>=2&contract==3) | (Q>=1&contract==2) | (Q>=1&contract==1) ) ; 
for ii=1:NActive+NMarginal1+NMarginal2 
    if ii<=NActive
        tempid = Activeids(ii) ; 
    elseif ii>NActive && ii<=NActive+NMarginal1
        tempid = Marginalids1(ii-NActive) ;
    else
        tempid = Marginalids2(ii-NActive-NMarginal1) ;
    end
    student_id_ltemp(ii) = tempid ;
    Qtemp(ii)            = Q(student_id_l==tempid) ;
    Ttemp(ii)            = T(student_id_l==tempid) ;
    ThetaEtemp(ii)       = ThetaE(student_id_l==tempid) ;
    contracttemp(ii)     = contract(student_id_l==tempid) ;
end
student_id_l = student_id_ltemp ;  %%%%This line now implies that student_id_l = [Activeids;Marginalids1;Marginalids2] 
Q            = Qtemp ;
T            = Ttemp ;
ThetaE       = ThetaEtemp ;
contract     = contracttemp ;
clear student_id_ltemp Qtemp Ttemp basetemp wagetemp contracttemp
if nargout==2
    out.student_id = student_id_l ;
    out.contract   = contract ;
    out.Q          = Q ;
    out.T          = T ;
    out.ThetaE     = ThetaE ;
end

Qbar    = length(SimTimeCum(:,1,1)) ; %%%%This is the maximum simulated output level (which may exceed the maximum observed value of 80)
Sims    = length(SimTimeCum(1,:,1)) ; %%%%This is the number of simulated work histories used for coarse search
if nargout==2
    Sims2 = Sims ;
    Sims3 = Sims ;
else
    Sims2   = max(10,min([Sims/4;40])) ; %%%%This is the number of simulated work histories used for fine search
    Sims3   = max(5,min([Sims2/4;10])) ; %%%%This is the number of simulated work histories used for counterfactual projections
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  BACKING OUT ThetaL, given ThetaE, Q, T:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%The idea behind the algorithm for backing out ThetaL is that we want to have mean simulated work time match T as closely as possible.


%%%%The following will be used as needed for extrapolation of the cost function in time space:
[Cend,Cprimeend,Cdblprimeend] = BsplineEval3(knotvec_c,pi_c,knotvec_c(end)) ; %%%%Upper tail levels and derivatives of the cost function


ThetaLACT     = nan(length(Activeids),2) ;  %%%%This will contain ThetaL values implied by the current parameter vector, for active students with Q<80
ThetaLMRG1    = nan(length(Marginalids1),2) ;  %%%%This will contain ThetaL values implied by the current parameter vector, for marginal students in contract 1
ThetaLMRG2    = nan(length(Marginalids2),2) ;  %%%%This will contain ThetaL values implied by the current parameter vector, for marginal students in contract 2
ThetaL80      = nan( sum(Q==80), length(Q80freq(:,1,1)) + 1 ) ;  %%%%This will contain a set of ThetaL values for each student with Q=80, based on extrapolated upper tail quantiles
ThetaL80(:,1) = student_id_l(Q==80) ;  %%%%This matrix will contain each student's id (only 
                                         %%%%students who achieved full output of Q=80) in the first 
                                         %%%%column, and the various possible values of ThetaL in 
                                         %%%%the subsequent columns, corresponding to the
                                         %%%%appropriate values of Q in Q80freq for that person.

%%%%The following are used for counterfactual computations
CFQsampACT_1_2 = nan(Sims3*sum(Q<80 & contract==1),1) ;  %%%%For simulating counterfactual choices of group-1 individuals under contract assingment 2
CFQsampACT_1_3 = CFQsampACT_1_2 ;  %%%%For simulating counterfactual choices of group-1 individuals under contract assingment 3
CFQsampACT_2_1 = nan(Sims3*sum(Q<80 & contract==2),1) ;  %%%%For simulating counterfactual choices of group-2 individuals under contract assingment 1
CFQsampACT_2_3 = CFQsampACT_2_1 ;  %%%%For simulating counterfactual choices of group-2 individuals under contract assingment 1
CFQsampACT_3_1 = nan(Sims3*sum(Q<80 & contract==3),1) ;  %%%%For simulating counterfactual choices of group-3 individuals under contract assingment 1
CFQsampACT_3_2 = CFQsampACT_3_1 ;  %%%%For simulating counterfactual choices of group-3 individuals under contract assingment 1
CFQsampMRG_1_2 = nan(Sims3*length(Marginalids1),1) ;  %%%%For simulating counterfactual choices of MARGINAL group-1 individuals under contract assingment 2
CFQsampMRG_1_3 = CFQsampMRG_1_2 ;  %%%%For simulating counterfactual choices of MARGINAL group-1 individuals under contract assingment 3
CFQsampMRG_2_1 = nan(Sims3*length(Marginalids2),1) ;  %%%%For simulating counterfactual choices of MARGINAL group-2 individuals under contract assingment 1
CFQsampMRG_2_3 = CFQsampMRG_2_1 ;  %%%%For simulating counterfactual choices of MARGINAL group-2 individuals under contract assingment 3

indxACT_1_2    = 0 ; indxACT_1_3 = 0 ;
indxACT_2_1    = 0 ; indxACT_2_3 = 0 ;
indxACT_3_1    = 0 ; indxACT_3_2 = 0 ;
indxMRG_1_2    = 0 ; indxMRG_1_3 = 0 ; 
indxMRG_2_1    = 0 ; indxMRG_2_3 = 0 ; 
%%%%Similar as above, but for folks with full output (i.e., 80 completed learning tasks)
CFQsamp80_1_2  = nan(Sims3*sum(Q==80 & contract==1)*sum(~isnan(Q80freq(:,1,1))),1) ;
CFQsamp80_1_3  = CFQsamp80_1_2 ;
CFQsamp80_2_1  = nan(Sims3*sum(Q==80 & contract==2)*sum(~isnan(Q80freq(:,1,2))),1) ;
CFQsamp80_2_3  = CFQsamp80_2_1 ;
CFQsamp80_3_1  = nan(Sims3*sum(Q==80 & contract==3)*sum(~isnan(Q80freq(:,1,3))),1) ;
CFQsamp80_3_2  = CFQsamp80_3_1 ;
indx80_1_2     = 0 ; indx80_1_3 = 0 ;
indx80_2_1     = 0 ; indx80_2_3 = 0 ;
indx80_3_1     = 0 ; indx80_3_2 = 0 ;




grosspay1     = pay(:,1)*ones(1,Sims) ;
grosspay2     = pay(:,2)*ones(1,Sims) ;
grosspay3     = pay(:,3)*ones(1,Sims) ;

ThetaLACTsave = nan(length(ThetaLACT),3) ;
for ii=1:length(Q)  %%%%Remember: on line 30 we set this up so that the first NActive values of
    %%%%Q/T/contract/etc. are for Activeids, and the next two blocks are for
    %%%%Marginalids1 and Marginalids2.  SimTimeCum was set up with the same
    %%%%sequencing for pages.
    if Q(ii)<80   %%%%This means the current student produced at least 1 unit, but less than full output
        tempcon  = contract(ii) ;  %%%%This is student ii's contract assignment
        if tempcon==1
            pi0 = 15 ;
            pi1 = 0.75 ;
        elseif tempcon==2
            pi0 = 10 ;
            pi1 = 1 ;
        else
            pi0 = 5 ;
            pi1 = 1.25 ;
        end
        temppage = find(student_id_l==student_id_l(ii)) ;  %%%IS THIS JUST THE SAME AS ii???
        %%%%This is the page in the array SimTimeCum that corresponds to the current student.
        %%%%Recall that student_id_l=[Activeids;Marginalids1;Marginalids2], and this indexing
        %%%%corresponds to the indexing of pages in SimTimeCum.
        tempQ      = Q(ii) ; %%%%This is student ii's OBSERVED cumulative output Q.
        tempT      = T(ii) ; %%%%This is student ii's OBSERVED cumulative work time T.
        Tdom2      = (0:3:max(120,tempT*1.5))' ;  %%%%This will be used for modelling continuation values.  We project 2 hours into future work, at 3 minute intervals...
        %%%%Here we retrieve the student's simulated work-time sequences.
        SimTime    = SimTimeCum(:,:,temppage) ; %#ok
        SimQuseful = max(21,min(tempQ*3,130)) ;  % %%%%For most people, 80+simulations is excessive, so we use this to see how much we need to match their ThetaL
        SimTime    = reshape(SimTime,Qbar*Sims,1) ;
        Cost       = nan(size(SimTime)) ;
        
        %%%%The following objects will be passed to the sub-routine OptStop.m:
        tempThetaE = ThetaE(ii) ; %%%%This is student ii's productivity type.
        tempseg    = find(cutoffs<tempThetaE,1,'last') ;
        if isempty(tempseg)
            tempseg = 1 ;
        elseif tempseg==length(cutoffs)
            tempseg = length(cutoffs) - 1 ;
        end
        
        %%%%The following conditional loop will be used as needed for extrapolation of the cost function in time space:
        if max(SimTime)<=knotvec_c(end) 
            Cost = BsplineEval3(knotvec_c,pi_c,SimTime) ;
        else  %%%%If we have to extrapolate, we use a 2nd-order taylor polynomial centered at knotvec_c(end)
            interpindx = SimTime <=knotvec_c(end) ;
            Cost(interpindx) =  BsplineEval3(knotvec_c,pi_c,SimTime(interpindx)) ;
            extrapindx = SimTime>knotvec_c(end) ;
            Cost(extrapindx) = Cend + Cprimeend*(SimTime(extrapindx)-knotvec_c(end)) + ...
                max(0,Cdblprimeend)*(SimTime(extrapindx)-knotvec_c(end)).^2 ;
        end
        SimTime     = reshape(SimTime,Qbar,Sims) ;
        Cost        = reshape(Cost,Qbar,Sims) ;
        
        SimTimetemp  = SimTime(1:SimQuseful,:) ;  %%%%The i-th row and j-th column of this contains the cumulative time to achieve i units of output, given the j-th sequence of simulated work-time shocks
        SimTimetemp2 = [zeros(1,Sims); SimTimetemp] ;
        Costtemp     = Cost(1:SimQuseful,:) ;  %%%%The i-th row and j-th column of this contains the cumulative cost (in levels) of achieving i units of output, given the j-th sequence of simulated work-time shocks
        
        %%%%The following two objects will be used to compute simulated profits:
        if tempcon==1
            grosspay = grosspay1(1:SimQuseful+1,:) ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        elseif tempcon==2
            grosspay = grosspay2(1:SimQuseful+1,:) ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        elseif tempcon==3
            grosspay = grosspay3(1:SimQuseful+1,:) ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        end
        commoncost   = [zeros(size(Costtemp(1,:))); Costtemp] ;    %%%%The dimension here is Qbar+1xSims, where the first row is gross profits/costs for doing nothing.        
        
        %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
        %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
        %%%%                           BISECTION SEARCH TO RECOVER THE Theta_L VALUES:
        %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
        %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
        
        %%%%%%%%11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111%%%%%%%%
        %%%%%%%%STEP 1, COARSE SEARCH: find a suitable search region from a broad set of possible values in ThetaL space.  To
        %%%%%%%%do this, we look for a value of ThetaLlo where the MYOPIC stopping point (a lower bound on the true stopping
        %%%%%%%%point) is the same as Ttemp, on average.  We also similarly fins a ThetaLhi value where the PERFECT FORESIGHT
        %%%%%%%%stopping point (an upper bound on the true stopping point) is the same as Ttemp, on average.  
        %%%%%%%%11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111%%%%%%%%  
        ThetaLMlo       = 10^-8 ; ThetaLMhi = 10 ; %%%%"M" is for MYOPIC.  related objects will be used to find the lower bound of the search region for step 2 (fine search)
        log10Mdiff      = log10(ThetaLMhi)-log10(ThetaLMlo) ; ThetaLMmid = 10^( log10(ThetaLMlo) + log10Mdiff/2 ) ; 
        ThetaL0Mgrid    = [ThetaLMlo; ThetaLMmid; ThetaLMhi] ; %%%%This contains both the end points and the midpoint
        tempTMtildemean = nan(size(ThetaL0Mgrid)) ; %%%%This will hold mean simulated myopic stop points in T-space for each element of ThetaL0grid.
        
        ThetaLPFlo       = 10^-8 ; ThetaLPFhi = 10 ; %%%%"PF" is for PERFECT FORESIGHT.  related objects will be used to find the lower bound of the search region for step 2 (fine search)
        log10PFdiff      = log10(ThetaLPFhi)-log10(ThetaLPFlo) ; ThetaLPFmid = 10^( log10(ThetaLPFlo) + log10PFdiff/2 ) ; 
        ThetaL0PFgrid    = [ThetaLPFlo; ThetaLPFmid; ThetaLPFhi] ; %%%%This contains both the end points and the midpoint 
        tempTPFtildemean = nan(size(ThetaL0PFgrid)) ;  %%%%This will hold mean simulated perfect-foresight stop points in T-space for each element of ThetaL0grid.

        iterM         = 0 ;
        iterPF        = 0 ;
        maxsearchiter = 25 ;
        searchtol1    = 10^-5 ;  %%%%This will be used for the search interval width
        searchtol2    = 10^-2 ;  %%%%This will be used for the objective function difference
        %%%%%%MYOPIC STOP POINTS TO REFINE THE LOWER BOUND OF THE COARSE SEARCH REGION
        for jj=2 %1:length(ThetaL0Mgrid)  %%%%This loop sets the initial values of the bisection search %%%%NOTE: SINCE WE ARE SURE ABOUT THE ENDPOINTS BEING WIDE ENOUGH, WE ONLY NEED THE OBJECTIVE FUNCTION VALUE AT THE MIDPOINT.
            profitmyopic = grosspay - ThetaL0Mgrid(jj)*commoncost ;
            profnegdiff = (profitmyopic(2:end,:) - profitmyopic(1:end-1,:)) < 0 ;  %%%%This is a matrix of 1s and 0s indicating whether successive total profits fell from one unit to the next.
                                                                                   %%%%We start from the second row, since myopic profits always drop from zero to one unit.
            [~,QstopMindx] = max(profnegdiff) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.  
            QstopMindx(QstopMindx==1 & sum(profnegdiff)==0) = length(profnegdiff(:,1));  %%%% EDGE CASE: When a given column of profnedgiff was all zeros, 
                                                                                           %%%% meaning profits trended positive throughtout the simulated shock sequence. In this 
                                                                                           %%%% case, QstopMindx(ss)==1 and profnegdiff(1,ss)==0. 
            tempTMtildemean(jj) = mean(SimTimetemp2(sub2ind(size(SimTimetemp2),QstopMindx,(1:Sims))))  ;
        end
        while log10Mdiff>searchtol1  && iterM<maxsearchiter  %%%%This loop takes the initial values and
            if abs(tempTMtildemean(2)-1.03*tempT)<=searchtol2  %%%%This is the secondary stopping criterion: if the midpoint T-value is very close to tempT, we take the midpoint, and it doesn't matter how wide the search region is...
                break %%%%This terminates the while loop
            end
            if tempTMtildemean(2) > 1.03*tempT  %%%%If the midpoint gives a T-value above tempT, we can throw away the lower endpoint, which would give a T-value even higher. 
                ThetaL0Mgrid(1)    = ThetaL0Mgrid(2) ;  %%%%Re-define the lower endpoint as the midpoint, and then re-define the midpoint below 
                tempTMtildemean(1) = tempTMtildemean(2) ;
            else  %%%%Otherwise, the midpoint gives a T-value below tempT, and we can throw away the upper endpoint, which would give a T-value even lower. 
                ThetaL0Mgrid(3)    = ThetaL0Mgrid(2) ;  %%%%Re-define the upper endpoint as the midpoint, and then re-define the midpoint below 
                tempTMtildemean(3) = tempTMtildemean(2) ;
            end
            log10Mdiff      = log10(ThetaL0Mgrid(3))-log10(ThetaL0Mgrid(1)) ;
            ThetaL0Mgrid(2) = 10^( log10(ThetaL0Mgrid(1)) + log10Mdiff/2 ) ;
            profitmyopic    = grosspay - ThetaL0Mgrid(2)*commoncost ;
            profnegdiff     = (profitmyopic(2:end,:) - profitmyopic(1:end-1,:)) < 0 ;  %%%%This is a matrix of 1s and 0s indicating whether successive total profits fell from one unit to the next.
            [~,QstopMindx]  = max(profnegdiff) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.
            QstopMindx(QstopMindx==1 & sum(profnegdiff)==0) = length(profnegdiff(:,1));  %%%% EDGE CASE: When a given column of profnedgiff was all zeros, 
                                                                                           %%%% meaning profits trended positive throughtout the simulated shock sequence. In this 
                                                                                           %%%% case, QstopMindx(ss)==1 and profnegdiff(1,ss)==0.  
            temp = sort(SimTimetemp2(sub2ind(size(SimTimetemp2),QstopMindx,(1:Sims)))) ;
            tempTMtildemean(2) = mean(temp) ;
            iterM = iterM + 1 ;
        end
        ThetaLlo = ThetaL0Mgrid(2) ;
        
        %%%%%%PERFECT FORESIGHT STOP POINTS TO REFINE THE UPPER BOUND OF THE COARSE SEARCH REGION
        for jj=2 %1:length(ThetaL0PFgrid)  %%%%This loop sets the initial values of the bisection search %%%%NOTE: SINCE WE ARE SURE ABOUT THE ENDPOINTS BEING WIDE ENOUGH, WE ONLY NEED THE OBJECTIVE FUNCTION VALUE AT THE MIDPOINT.
            profitmyopic = grosspay - ThetaL0PFgrid(jj)*commoncost ;
            [~,QstopPFindx] = max(profitmyopic) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.  
            tempTPFtildemean(jj) = mean(SimTimetemp2(sub2ind(size(SimTimetemp2),QstopPFindx,(1:Sims)))) ; 
        end
        while log10PFdiff>searchtol1  && iterPF<maxsearchiter  %%%%This loop takes the initial values and
            if abs(tempTPFtildemean(2)-tempT)<=searchtol2  %%%%This is the secondary stopping criterion: if the midpoint T-value is very close to tempT, we take the midpoint, and it doesn't matter how wide the search region is...
                break %%%%This terminates the while loop
            end
            if tempTPFtildemean(2) > tempT  %%%%If the midpoint gives a T-value above tempT, we can throw away the lower endpoint, which would give a T-value even higher. 
                ThetaL0PFgrid(1)    = ThetaL0PFgrid(2) ;  %%%%Re-define the lower endpoint as the midpoint, and then re-define the midpoint below 
                tempTPFtildemean(1) = tempTPFtildemean(2) ;
            else  %%%%Otherwise, the midpoint gives a T-value below tempT, and we can throw away the upper endpoint, which would give a T-value even lower. 
                ThetaL0PFgrid(3)    = ThetaL0PFgrid(2) ;  %%%%Re-define the upper endpoint as the midpoint, and then re-define the midpoint below 
                tempTPFtildemean(3) = tempTPFtildemean(2) ;
            end
            log10PFdiff      = log10(ThetaL0PFgrid(3))-log10(ThetaL0PFgrid(1)) ;
            ThetaL0PFgrid(2) = 10^( log10(ThetaL0PFgrid(1)) + log10PFdiff/2 ) ;
            profitmyopic     = grosspay - ThetaL0PFgrid(jj)*commoncost ;
            [~,QstopPFindx]  = max(profitmyopic) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined. 
            temp = sort(SimTimetemp2(sub2ind(size(SimTimetemp2),QstopPFindx,(1:Sims)))) ;
            tempTPFtildemean(2) = mean(temp) ;
            iterPF = iterPF + 1 ;
        end
        ThetaLhi = ThetaL0PFgrid(2) ;
        %%%%%%%%22222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222%%%%%%%%
        %%%%%%%%STEP 2, FINE SEARCH: using the values of ThetaLlo and ThetaLhi above, we finish the inversion routine to 
        %%%%%%%%recover ThetaL by evaluating the true mean stop time at the endpoints and two intermediate points, and then
        %%%%%%%%interpolating tempT. In doing so, we economize on continuation function evaluations (the slow part) by solving for
        %%%%%%%%myopic stop times first, and then taking incremental steps forward from the stopping point using the full
        %%%%%%%%value function, including continuation value.  We stop once the value function difference goes negative from
        %%%%%%%%one unit to the next, as this corresponds to the stopping point (in Q space) of the full dynamic program.
        %%%%%%%%22222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222%%%%%%%%
        %%%%IN ORDER TO ECONOMIZE ON COMPUTATION TIME AND AVOID TOO MANY VALUE FUNCTION EVALUATIONS, WE PERFORM FINE SEARCH
        %%%%ON A STRICT SUBSET OF THE SHOCK SEQUENCES.  BUT, IN ORDER TO GET MAXIMAL STATISTICAL POWER, WE CHOOSE A 
        %%%%REPRESENTATIVE SUBSET OF SEQUENCES AS THE (1/Sims2) QUANTILES OF SimTime(upshift*tempQ,:), WHERE upshift IS 
        %%%%DEFINED BELOW:  
        if nargout==2
            Sims2indx = (1:Sims2) ;
        else
            upshift = 1+std(SimTimetemp(tempQ,:))/mean(SimTimetemp(tempQ,:)) ;
            SimTimes_attempQ      = sort(SimTimetemp(round(upshift*tempQ),:)) ;
            SimTimes_attempQshort = SimTimes_attempQ(Sims/(2*Sims2):Sims/Sims2:Sims-(Sims/(2*Sims2))) ; %%%%This line picks out the (1/Sims2) quantiles of the shock sequences.
            Sims2indx             = nan(1,Sims2) ;
            for kk=1:Sims2
                Sims2indx(kk) = find(SimTimetemp(round(upshift*tempQ),:)==SimTimes_attempQshort(kk)) ;
            end
            Sims2indx = sort(Sims2indx) ;
        end

        ThetaL0grid = [ThetaLlo; 10^( (log10(ThetaLlo)+log10(ThetaLhi))/2 ); ThetaLhi] ;
        tempTtildemean = nan(size(ThetaL0grid)) ; 
        for jj=1:length(ThetaL0grid)  %%%%This loop sets the initial values of the bisection search %%%%NOTE: SINCE WE ARE SURE ABOUT THE ENDPOINTS BEING WIDE ENOUGH, WE ONLY NEED THE OBJECTIVE FUNCTION VALUE AT THE MIDPOINT.
            profitmyopic = grosspay(:,Sims2indx) - ThetaL0grid(jj)*commoncost(:,Sims2indx) ;
            profnegdiff  = (profitmyopic(2:end,:) - profitmyopic(1:end-1,:)) < 0 ; %%%%This is a matrix of 1s and 0s indicating whether successive total profits fell from one unit to the next.
                                                                                   %%%%We start from the second row, since myopic profits always drop from zero to one unit.
            [~,Qstopindx] = max(profnegdiff) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.  
            Qstopindx(Qstopindx==1 & sum(profnegdiff)==0) = length(profnegdiff(:,1));  %%%% EDGE CASE: When a given column of profnedgiff was all zeros, 
                                                                                           %%%% meaning profits trended positive throughtout the simulated shock sequence. In this 
                                                                                           %%%% case, QstopMindx(ss)==1 and profnegdiff(1,ss)==0, and we re-set the value to Qbar.     
            ThetaEweight = (log(tempThetaE)-log(cutoffs(tempseg+1)))/(log(cutoffs(tempseg))-log(cutoffs(tempseg+1))) ;  %%%%This is the weight that will be placed on the lower-bound (I.E, MOST PRODUCTIVE TYPE) of the ThetaE segment.
            EQgivenT = [zeros(1,Qbar+1); exp( ThetaEweight*log(EQgivenT_ub(2:end,:,tempseg))+(1-ThetaEweight)*log(EQgivenT_lb(2:end,:,tempseg)) )] ;
            tempTtilde = zeros(size(profitmyopic(1,:))) ;

            for ss=1:Sims2 
                Qcurrent = Qstopindx(ss)-1 ;
                if Qcurrent<SimQuseful  %%%%If the current simulated stop time is already the upper bound, we just take the simulated upper bound work time.
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%INITIATE THE VALUE FUNCTION AT THE CURRENT Q:
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    if Qcurrent==0
                        Tdom3 = Tdom2 ;  %%%%THIS WILL BE USED LATER TO COMPUTE COSTS, STARTING AT THE CURRENT TIME COMMITMENT BASELINE.
                    else
                        Tdom3 = SimTimetemp(Qcurrent,Sims2indx(ss)) + Tdom2 ;  %%%%THIS WILL BE USED LATER TO COMPUTE COSTS.
                    end
                    %%%THIS IS THE CONTINUATION REVENUE OF FUTURE WORK FROM THE CURRENT Q VALUE:
                    CRev = interp1(Tdom,EQgivenT(:,Qcurrent+1),Tdom2,'linear') ; %%%%START WITH EXPECTED OUTPUT q FOR EACH VALUE OF t IN Tdom3
                    if Qcurrent==0
                        CRev(CRev<2)  = (CRev(CRev<2)/2)*(2*pi1 + pi0) ; 
                        CRev(CRev>=2) = CRev(CRev>=2)*pi1 + pi0 ; 
                    elseif Qcurrent==1
                        CRev(CRev<1)  = CRev(CRev<1)*(2*pi1 + pi0) ; 
                        CRev(CRev>=1) = CRev(CRev>=1)*pi1 + pi1 + pi0 ; 
                    else
                        CRev = CRev*pi1 ; 
                    end
                    if max(Tdom3)<=knotvec_c(end)  %%%%EDGE CASE: Need some special adjustments if Tdom3 goes outside the usual cost domain
                        CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qcurrent+1,Sims2indx(ss)) ; %%%THIS IS THE CONTINUATION COST OF FUTURE WORK  NOTE: commoncost has a first row of zeros for zero output
                    else
                        CCost                   = zeros(size(Tdom3)) ;
                        interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                        CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                        CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +... %%%%For extrapolation, we use a 2nd-order Taylor polynomial
                                                         (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                        CCost                   = CCost - commoncost(Qcurrent+1,Sims2indx(ss)) ; %%%NOTE: commoncost has a first row of zeros for zero output
                    end
                    CValcurrent   = max(0,max(CRev - ThetaL0grid(jj)*CCost)) ;
                    ValFuncurrent = profitmyopic(Qcurrent+1,ss) + CValcurrent ;
                    stopflag = 0 ;  %%%%This will be used to signal when we can stop stepping forward looking for the stop point.
                    while stopflag==0
                        if CValcurrent==0
                            break
                        end
                        Qnext = Qcurrent + 1 ;
                        Tdom3 = SimTimetemp(Qnext,Sims2indx(ss)) + Tdom2 ;
                        %%%THIS IS THE CONTINUATION REVENUE OF FUTURE WORK:
                        CRev = interp1(Tdom,EQgivenT(:,Qnext+1),Tdom2,'linear') ; %%%%START WITH EXPECTED OUTPUT q FOR EACH VALUE OF t IN Tdom3
                        if Qnext==1
                            CRev(CRev<1)  = CRev(CRev<1)*(2*pi1 + pi0) ;
                            CRev(CRev>=1) = CRev(CRev>=1)*pi1 + pi1 + pi0 ;
                        else
                            CRev = CRev*pi1 ;
                        end
                        if max(Tdom3)<=knotvec_c(end)  %%%%EDGE CASE: Need some special adjustments if Tdom3 goes outside the usual cost domain
                            CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qnext+1,Sims2indx(ss)) ; %%%THIS IS THE CONTINUATION COST OF FUTURE WORK  NOTE: commoncost has a first row of zeros for zero output
                        else
                            CCost                   = zeros(size(Tdom3)) ;
                            interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                            CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                            CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +... %%%%For extrapolation, we use a 2nd-order Taylor polynomial
                                                             (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                            CCost                   = CCost  - commoncost(Qcurrent+1,Sims2indx(ss)) ; %%%NOTE: commoncost has a first row of zeros for zero output
                        end
                        CValnext   = max(0,max(CRev - ThetaL0grid(jj)*CCost)) ;
                        ValFunnext = profitmyopic(Qnext+1,ss) + CValnext ;
                        if  ValFuncurrent > ValFunnext  %%%%This means we stop now... Also, if the Continuation Value is zero, there is no need to go any further, because myopic profits and the value function are one and the same.
                            stopflag = 1 ;
                        else
                            Qcurrent      = Qnext ;
                            CValcurrent   = CValnext ;
                            ValFuncurrent = ValFunnext ;
                        end
                        if Qcurrent+1>=SimQuseful  %%%%If the next step would take us to/above the simulated work upper bound, we just stop.
                            stopflag = 1 ;
                        end

                    end

                    if Qcurrent==0
                        tempTtilde(ss) = 0 ;
                    else
                        tempTtilde(ss) = SimTimetemp(Qcurrent,Sims2indx(ss)) ;
                    end
                else
                    tempTtilde(ss) = SimTimetemp(Qcurrent,Sims2indx(ss)) ;
                end

            end
            tempTtildemean(jj) = mean(tempTtilde)  ; 
        end  

        if tempQ>=2 
            if tempT>tempTtildemean(1)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL    = ThetaL0grid(1) ; 
                ThetaLACT(ii,1) = tempThetaL ; ThetaLACT(ii,2) = student_id_l(ii) ; 
            elseif tempT<tempTtildemean(end)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL    = ThetaL0grid(end) ; 
                ThetaLACT(ii,1) = tempThetaL ; ThetaLACT(ii,2) = student_id_l(ii) ; 
            else  %%%% Otherwise, we interpolate, taking into account a final edge case were multiple simulated means may be the same...
                if tempTtildemean(1)==tempTtildemean(2) && tempTtildemean(2)==tempTtildemean(3)
                    tempThetaL = ThetaL0grid(2) ;
                    ThetaLACT(ii,1) =tempThetaL ;  ThetaLACT(ii,2) = student_id_l(ii) ; 
                elseif tempTtildemean(1)==tempTtildemean(2) 
                    tempThetaL = 10^(interp1(tempTtildemean(2:3),log10(ThetaL0grid(2:3)),tempT,'linear')) ;
                    ThetaLACT(ii,1) = tempThetaL ; ThetaLACT(ii,2) = student_id_l(ii) ; 
                elseif tempTtildemean(2)==tempTtildemean(3) 
                    tempThetaL = 10^(interp1(tempTtildemean(1:2),log10(ThetaL0grid(1:2)),tempT,'linear')) ;
                    ThetaLACT(ii,1) = tempThetaL ; ThetaLACT(ii,2) = student_id_l(ii) ; 
                else
                    tempThetaL = 10^(interp1(tempTtildemean,log10(ThetaL0grid),tempT,'linear')) ;
                    ThetaLACT(ii,1) = tempThetaL ; ThetaLACT(ii,2) = student_id_l(ii) ; 
                end
            end 
            ThetaLACTsave(ii,:) = ThetaL0grid ;
        elseif tempQ==1 && tempcon==1
            if tempT>tempTtildemean(1)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL = ThetaL0grid(1) ;
                ThetaLMRG1(ii-NActive,1) =  tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
            elseif tempT<tempTtildemean(end)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL = ThetaL0grid(end) ;
                ThetaLMRG1(ii-NActive,1) = tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
            else  %%%% Otherwise, we interpolate, taking into account a final edge case were multiple simulated means may be the same...
                if tempTtildemean(1)==tempTtildemean(2) && tempTtildemean(2)==tempTtildemean(3)
                    tempThetaL = ThetaL0grid(2) ;
                    ThetaLMRG1(ii-NActive,1) = tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
                elseif tempTtildemean(1)==tempTtildemean(2) 
                    tempThetaL = 10^(interp1(tempTtildemean(2:3),log10(ThetaL0grid(2:3)),tempT,'linear')) ;
                    ThetaLMRG1(ii-NActive,1) = tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
                elseif tempTtildemean(2)==tempTtildemean(3) 
                    tempThetaL = 10^(interp1(tempTtildemean(1:2),log10(ThetaL0grid(1:2)),tempT,'linear')) ;
                    ThetaLMRG1(ii-NActive,1) = tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
                else
                    tempThetaL = 10^(interp1(tempTtildemean,log10(ThetaL0grid),tempT,'linear')) ;
                    ThetaLMRG1(ii-NActive,1) = tempThetaL ; ThetaLMRG1(ii-NActive,2) = student_id_l(ii) ; 
                end
            end
        elseif tempQ==1 && tempcon==2
            if tempT>tempTtildemean(1)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL = ThetaL0grid(1) ; 
                ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
            elseif tempT<tempTtildemean(end)  %%%%This means that the observed action is outside the simulated range.  Since coarse search is based on a larger sample, we use those values as bounds.
                tempThetaL = ThetaL0grid(end) ;
                ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
            else  %%%% Otherwise, we interpolate, taking into account a final edge case were multiple simulated means may be the same...
                if tempTtildemean(1)==tempTtildemean(2) && tempTtildemean(2)==tempTtildemean(3)
                    tempThetaL =ThetaL0grid(2) ;
                    ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
                elseif tempTtildemean(1)==tempTtildemean(2) 
                    tempThetaL = 10^(interp1(tempTtildemean(2:3),log10(ThetaL0grid(2:3)),tempT,'linear')) ;
                    ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
                elseif tempTtildemean(2)==tempTtildemean(3) 
                    tempThetaL = 10^(interp1(tempTtildemean(1:2),log10(ThetaL0grid(1:2)),tempT,'linear')) ;
                    ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
                else 
                    tempThetaL = 10^(interp1(tempTtildemean,log10(ThetaL0grid),tempT,'linear')) ;
                    ThetaLMRG2(ii-NActive-NMarginal1,1) = tempThetaL ; ThetaLMRG2(ii-NActive-NMarginal1,2) = student_id_l(ii) ; 
                end
            end
        end
        %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
        %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
        %%%%%%%%                           STEP 3, COUNTERFACTUAL SIMULATIONS
        %%%%%%%%  NOTE: for descriptive code comments, see Step 2 above.
        %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
        %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
        %%%%IN ORDER TO ECONOMIZE ON COMPUTATION TIME AND AVOID TOO MANY VALUE FUNCTION EVALUATIONS, WE PERFORM FINE SEARCH
        %%%%ON A STRICT SUBSET OF THE SHOCK SEQUENCES.  BUT, IN ORDER TO GET MAXIMAL STATISTICAL POWER, WE CHOOSE A
        %%%%REPRESENTATIVE SUBSET OF SEQUENCES AS THE (1/Sims2) QUANTILES OF SimTime(upshift*tempQ,:), WHERE upshift IS
        %%%%DEFINED BELOW:
        if nargout==2
            Sims3indx = (1:Sims3) ;
            commoncost   = [zeros(size(Cost(1,:))); Cost] ;    %%%%The dimension here is Qbar+1xSims, where the first row is gross profits/costs for doing nothing.
        else
            SimTimes_attempQshort = SimTimes_attempQ(Sims/(2*Sims3):Sims/Sims3:Sims-(Sims/(2*Sims3))) ; %%%%This line picks out the (1/Sims2) quantiles of the shock sequences.
            Sims3indx             = nan(1,Sims3) ;
            commoncost   = [zeros(size(Cost(1,:))); Cost] ;    %%%%The dimension here is Qbar+1xSims, where the first row is gross profits/costs for doing nothing.
            for kk=1:Sims3
                Sims3indx(kk) = find(SimTime(round(upshift*tempQ),:)==SimTimes_attempQshort(kk)) ;
            end
            Sims3indx = sort(Sims3indx) ;
        end
        if tempcon==1
            profitmyopicCF1 = grosspay2(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            profitmyopicCF2 = grosspay3(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
        elseif tempcon==2
            profitmyopicCF1 = grosspay1(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            profitmyopicCF2 = grosspay3(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
        elseif tempcon==3
            profitmyopicCF1 = grosspay1(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            profitmyopicCF2 = grosspay2(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
        end

        profnegdiffCF1  = (profitmyopicCF1(2:end,:) - profitmyopicCF1(1:end-1,:)) < 0 ; 
        [~,QstopindxCF1] = max(profnegdiffCF1) ; 
        QstopindxCF1(QstopindxCF1==1 & sum(profnegdiffCF1)==0) = length(profnegdiffCF1(:,1)); 
        tempQtildeCF1 = zeros(size(profitmyopicCF1(1,:))) ;

        profnegdiffCF2  = (profitmyopicCF2(2:end,:) - profitmyopicCF2(1:end-1,:)) < 0 ; 
        [~,QstopindxCF2] = max(profnegdiffCF2) ; 
        QstopindxCF2(QstopindxCF2==1 & sum(profnegdiffCF2)==0) = length(profnegdiffCF2(:,1));  
        tempQtildeCF2 = zeros(size(profitmyopicCF2(1,:))) ;
        for ss=1:Sims3 
            for kk=1:2  %%%%The value of kk indicates whether we are doing CF1 (the low piece-rate CF) or CF2 (the high piece-rate CF)
                if kk==1  %%%%This means we are currently doing the low-piece-rate CF
                    Qcurrent = QstopindxCF1(ss)-1 ;
                    if tempcon==1  %%%%In this case, we are doing the low-piece-rate CF which is contract 2, since the actual contract is 1
                        pi0CF = 10 ;
                        pi1CF = 1 ;
                    elseif tempcon>1 %%%%In this case, we are doing the low-piece-rate CF which is contract 1, since the actual contract is either 2 or 3
                        pi0CF = 15 ;
                        pi1CF = 0.75 ;
                    end
                elseif kk==2  %%%%This means we are currently doing the high-piece-rate CF
                    Qcurrent = QstopindxCF2(ss)-1 ;
                    if tempcon==3  %%%%In this case, we are doing the high-piece-rate CF which is contract 2, since the actual contract is 3
                        pi0CF = 10 ;
                        pi1CF = 1 ;
                    elseif tempcon<3 %%%%In this case, we are doing the high-piece-rate CF which is contract 3, since the actual contract is either 1 or 2
                        pi0CF = 5 ;
                        pi1CF = 1.25 ;
                    end
                end
                if Qcurrent<81
                    if Qcurrent==0
                        Tdom3 = Tdom2 ;
                    else
                        Tdom3 = SimTime(Qcurrent,Sims3indx(ss)) + Tdom2 ;
                    end
                    CRev = interp1(Tdom,EQgivenT(:,Qcurrent+1),Tdom2,'linear') ;
                    if Qcurrent==0
                        CRev(CRev<2)  = (CRev(CRev<2)/2)*(2*pi1CF + pi0CF) ;
                        CRev(CRev>=2) = CRev(CRev>=2)*pi1CF + pi0CF ;
                    elseif Qcurrent==1
                        CRev(CRev<1)  = CRev(CRev<1)*(2*pi1CF + pi0CF) ;
                        CRev(CRev>=1) = CRev(CRev>=1)*pi1CF + pi1CF + pi0CF ;
                    else
                        CRev = CRev*pi1CF ;
                    end
                    if max(Tdom3)<=knotvec_c(end)
                        CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qcurrent+1,Sims3indx(ss)) ;
                    else
                        CCost                   = zeros(size(Tdom3)) ;
                        interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                        CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                        CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +... 
                                                         (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                        CCost                   = CCost - commoncost(Qcurrent+1,Sims3indx(ss)) ; 
                    end
                    CValcurrent   = max(0,max(CRev - tempThetaL*CCost)) ;
                    if kk==1
                        ValFuncurrent = profitmyopicCF1(Qcurrent+1,ss) + CValcurrent ;
                    elseif kk==2
                        ValFuncurrent = profitmyopicCF2(Qcurrent+1,ss) + CValcurrent ;
                    end
                    stopflag = 0 ;  
                    while stopflag==0
                        if CValcurrent==0
                            break
                        end
                        Qnext = Qcurrent + 1 ;
                        Tdom3 = SimTime(Qnext,Sims3indx(ss)) + Tdom2 ;
                        CRev  = interp1(Tdom,EQgivenT(:,Qnext+1),Tdom2,'linear') ; %%%%START WITH EXPECTED OUTPUT q FOR EACH VALUE OF t IN Tdom3
                        if Qnext==1
                            CRev(CRev<1)  = CRev(CRev<1)*(2*pi1CF + pi0CF) ;
                            CRev(CRev>=1) = CRev(CRev>=1)*pi1CF + pi1CF + pi0CF ;
                        else
                            CRev = CRev*pi1CF ;
                        end
                        if max(Tdom3)<=knotvec_c(end)  
                            CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qnext+1,Sims3indx(ss)) ; 
                        else
                            CCost                   = zeros(size(Tdom3)) ;
                            interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                            CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                            CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +... 
                                                             (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                            CCost                   = CCost  - commoncost(Qcurrent+1,Sims3indx(ss)) ; 
                        end
                        CValnext   = max(0,max(CRev - tempThetaL*CCost)) ;
                        if kk==1
                            ValFunnext = profitmyopicCF1(Qnext+1,ss) + CValnext ;
                        elseif kk==2
                            ValFunnext = profitmyopicCF2(Qnext+1,ss) + CValnext ;
                        end
                        if  ValFuncurrent > ValFunnext  
                            stopflag = 1 ;
                        else
                            Qcurrent      = Qnext ;
                            CValcurrent   = CValnext ;
                            ValFuncurrent = ValFunnext ;
                        end
                        if Qcurrent+1>=81  
                            stopflag = 1 ;
                        end
                    end

                    if kk==1
                        tempQtildeCF1(ss) = Qcurrent ;
                    elseif kk==2
                        tempQtildeCF2(ss) = Qcurrent ;
                    end
                else
                    if kk==1
                        tempQtildeCF1(ss) = Qcurrent ;
                    elseif kk==2
                        tempQtildeCF2(ss) = Qcurrent ;
                    end
                end
            end
        end

        if tempQ>=2
            if tempcon==1
                CFQsampACT_1_2(indxACT_1_2+1:indxACT_1_2+Sims3) = tempQtildeCF1 ; indxACT_1_2 = indxACT_1_2 + Sims3 ;
                CFQsampACT_1_3(indxACT_1_3+1:indxACT_1_3+Sims3) = tempQtildeCF2 ; indxACT_1_3 = indxACT_1_3 + Sims3 ;
            elseif tempcon==2
                CFQsampACT_2_1(indxACT_2_1+1:indxACT_2_1+Sims3) = tempQtildeCF1 ; indxACT_2_1 = indxACT_2_1 + Sims3 ;
                CFQsampACT_2_3(indxACT_2_3+1:indxACT_2_3+Sims3) = tempQtildeCF2 ; indxACT_2_3 = indxACT_2_3 + Sims3 ;
            elseif tempcon==3
                CFQsampACT_3_1(indxACT_3_1+1:indxACT_3_1+Sims3) = tempQtildeCF1 ; indxACT_3_1 = indxACT_3_1 + Sims3 ;
                CFQsampACT_3_2(indxACT_3_2+1:indxACT_3_2+Sims3) = tempQtildeCF2 ; indxACT_3_2 = indxACT_3_2 + Sims3 ;
            end
        elseif tempQ==1
            if tempcon==1
                CFQsampMRG_1_2(indxMRG_1_2+1:indxMRG_1_2+Sims3) = tempQtildeCF1 ; indxMRG_1_2 = indxMRG_1_2 + Sims3 ;
                CFQsampMRG_1_3(indxMRG_1_3+1:indxMRG_1_3+Sims3) = tempQtildeCF2 ; indxMRG_1_3 = indxMRG_1_3 + Sims3 ;
            elseif tempcon==2
                CFQsampMRG_2_1(indxMRG_2_1+1:indxMRG_2_1+Sims3) = tempQtildeCF1 ; indxMRG_2_1 = indxMRG_2_1 + Sims3 ;
                CFQsampMRG_2_3(indxMRG_2_3+1:indxMRG_2_3+Sims3) = tempQtildeCF2 ; indxMRG_2_3 = indxMRG_2_3 + Sims3 ;
            end
        end
    end


end 
logThetaLlowt = mean((log10(ThetaLACT(Q>=50))-log10(ThetaLACTsave(Q>=50,3)))./...
                     (log10(ThetaLACTsave(Q>=50,1))-log10(ThetaLACTsave(Q>=50,3))),'omitnan') ;  %%%%We chose above 50 because this is roughly one standard deviation above the mean, conditional on Active status.

for ii=1:length(Q)  %%%%Remember: on line 30 we set this up so that the first NActive values of
    %%%%Q/T/contract/etc. are for Activeids, and the next two blocks are for
    %%%%Marginalids1 and Marginalids2.  SimTimeCum was set up with the same
    %%%%sequencing for pages.
    if Q(ii)==80
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%This part of the loop covers people in the upper mass point at 80 units of output.
        %%%%For these folks, to back out ThetaL values, we just match on extrapolated values of Q, conditional on their
        %%%%ThetaE type.  This is because we don't have extrapolated T values for them; only extrapolated Q values...
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        tempcon  = contract(ii) ;  %%%%This is student ii's contract assignment

        temppage = find(student_id_l==student_id_l(ii)) ;  %%%IS THIS JUST THE SAME AS ii???
        %%%%This is the page in the array SimTimeCum that corresponds to the current student.
        %%%%Recall that student_id_l=[Activeids;Marginalids1;Marginalids2], and this indexing
        %%%%corresponds to the indexing of pages in SimTimeCum.
        tempQ      = Q(ii) ; %%%%This is student ii's OBSERVED cumulative output Q.
        tempT      = T(ii) ; %%%%This is student ii's OBSERVED cumulative work time T.
        Tdom2      = (0:3:max(120,tempT*1.5))' ;  %%%%This will be used for modelling continuation values.  We project 2 hours into future work, at 3 minute intervals...
        %%%%Here we retrieve the student's simulated work-time sequences.
        SimTime    = SimTimeCum(:,:,temppage) ; %#ok
        SimTime    = reshape(SimTime,Qbar*Sims,1) ;
        Cost       = nan(size(SimTime)) ;
        

        
        if max(SimTime)<=knotvec_c(end)
            Cost = BsplineEval3(knotvec_c,pi_c,SimTime) ;
        else  %%%%If we have to extrapolate, we use a 2nd-order taylor polynomial centered at knotvec_c(end)
            interpindx = SimTime <=knotvec_c(end) ;
            Cost(interpindx) =  BsplineEval3(knotvec_c,pi_c,SimTime(interpindx)) ;
            extrapindx = SimTime>knotvec_c(end) ;
            Cost(extrapindx) = Cend + Cprimeend*(SimTime(extrapindx)-knotvec_c(end)) + ...
                max(0,Cdblprimeend)*(SimTime(extrapindx)-knotvec_c(end)).^2 ;
        end
        SimTime     = reshape(SimTime,Qbar,Sims) ;
        Cost        = reshape(Cost,Qbar,Sims) ;
        Q80temp = Q80freq(:,1,tempcon) ; Q80temp(isnan(Q80temp)) = [] ;  %%%%This retrieves the frequency table of extrapolated upper-tail quantiles

        
        %%%%The following two objects will be used to compute simulated profits:
        if tempcon==1
            grosspay = grosspay1 ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        elseif tempcon==2
            grosspay = grosspay2 ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        elseif tempcon==3
            grosspay = grosspay3 ;%%%%Remember that the pay matrix starts (first row) at pay for zero work...
        end
        commoncost   = [zeros(size(Cost(1,:))); Cost] ;    %%%%The dimension here is Qbar+1xSims, where the first row is gross profits/costs for doing nothing.
        ThetaLlo = nan(length(Q80temp),1) ;
        ThetaLhi = nan(length(Q80temp),1) ;
        %%%%The following four lines will be for counterfactual simulations...
        %%%%IN ORDER TO ECONOMIZE ON COMPUTATION TIME AND AVOID TOO MANY VALUE FUNCTION EVALUATIONS, WE PERFORM FINE SEARCH
        %%%%ON A STRICT SUBSET OF THE SHOCK SEQUENCES.  BUT, IN ORDER TO GET MAXIMAL STATISTICAL POWER, WE CHOOSE A
        %%%%REPRESENTATIVE SUBSET OF SEQUENCES AS THE (1/Sims2) QUANTILES OF SimTime(upshift*tempQ,:), WHERE upshift IS
        %%%%DEFINED BELOW:
        if nargout==2
            Sims3indx = (1:Sims3) ;
        else
            upshift = 1+std(SimTime(tempQ,:))/mean(SimTime(tempQ,:)) ;
            SimTimes_attempQ      = sort(SimTime(round(upshift*tempQ),:)) ;
            SimTimes_attempQshort = SimTimes_attempQ(Sims/(2*Sims3):Sims/Sims3:Sims-(Sims/(2*Sims3))) ; %%%%This line picks out the (1/Sims2) quantiles of the shock sequences.
            Sims3indx             = nan(1,Sims3) ;
            for kk=1:Sims3
                Sims3indx(kk) = find(SimTime(round(upshift*tempQ),:)==SimTimes_attempQshort(kk)) ;
            end
            Sims3indx = sort(Sims3indx) ;
        end
        for ll=1:length(Q80temp)
            tempQ = Q80temp(ll) ;
            %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
            %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
            %%%%                           BISECTION SEARCH TO RECOVER THE Theta_L VALUES:
            %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%
            %%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%BISECTION%%%%%%%%%%%

            %%%%%%%%11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111%%%%%%%%
            %%%%%%%%STEP 1, COARSE SEARCH: find a suitable search region from a broad set of possible values in ThetaL space.  To
            %%%%%%%%do this, we look for a value of ThetaLlo where the MYOPIC stopping point (a lower bound on the true stopping
            %%%%%%%%point) is the same as Ttemp, on average.  We also similarly fins a ThetaLhi value where the PERFECT FORESIGHT
            %%%%%%%%stopping point (an upper bound on the true stopping point) is the same as Ttemp, on average.
            %%%%%%%%11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111%%%%%%%%
            ThetaLMlo       = 10^-8 ; ThetaLMhi = 10 ; %%%%"M" is for MYOPIC.  related objects will be used to find the lower bound of the search region for step 2 (fine search)
            log10Mdiff      = log10(ThetaLMhi)-log10(ThetaLMlo) ; ThetaLMmid = 10^( log10(ThetaLMlo) + log10Mdiff/2 ) ;
            ThetaL0Mgrid    = [ThetaLMlo; ThetaLMmid; ThetaLMhi] ; %%%%This contains both the end points and the midpoint
            tempQMtildemean = nan(size(ThetaL0Mgrid)) ; %%%%This will hold mean simulated myopic stop points in T-space for each element of ThetaL0grid.

            ThetaLPFlo       = 10^-8 ; ThetaLPFhi = 10 ; %%%%"PF" is for PERFECT FORESIGHT.  related objects will be used to find the lower bound of the search region for step 2 (fine search)
            log10PFdiff      = log10(ThetaLPFhi)-log10(ThetaLPFlo) ; ThetaLPFmid = 10^( log10(ThetaLPFlo) + log10PFdiff/2 ) ;
            ThetaL0PFgrid    = [ThetaLPFlo; ThetaLPFmid; ThetaLPFhi] ; %%%%This contains both the end points and the midpoint
            tempQPFtildemean = nan(size(ThetaL0PFgrid)) ;  %%%%This will hold mean simulated perfect-foresight stop points in T-space for each element of ThetaL0grid.

            iterM         = 0 ;
            iterPF        = 0 ;
            maxsearchiter = 25 ;
            searchtol1    = 10^-5 ;  %%%%This will be used for the search interval width
            searchtol2    = 10^-2 ;  %%%%This will be used for the objective function difference
            %%%%%%MYOPIC STOP POINTS TO REFINE THE LOWER BOUND OF THE COARSE SEARCH REGION
            for jj=2 %1:length(ThetaL0Mgrid)  %%%%This loop sets the initial values of the bisection search %%%%NOTE: SINCE WE ARE SURE ABOUT THE ENDPOINTS BEING WIDE ENOUGH, WE ONLY NEED THE OBJECTIVE FUNCTION VALUE AT THE MIDPOINT.
                profitmyopic = grosspay - ThetaL0Mgrid(jj)*commoncost ;
                profnegdiff = (profitmyopic(2:end,:) - profitmyopic(1:end-1,:)) < 0 ;  %%%%This is a matrix of 1s and 0s indicating whether successive total profits fell from one unit to the next.
                                                                                       %%%%We start from the second row, since myopic profits always drop from zero to one unit.
                [~,QstopMindx] = max(profnegdiff) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.
                QstopMindx(QstopMindx==1 & sum(profnegdiff)==0) = length(profnegdiff(:,1));  %%%% EDGE CASE: When a given column of profnedgiff was all zeros,
                                                                                             %%%% meaning profits trended positive throughtout the simulated shock sequence. In this
                                                                                             %%%% case, QstopMindx(ss)==1 and profnegdiff(1,ss)==0.
                tempQMtildemean(jj) = mean(QstopMindx-1)  ;
            end
            while log10Mdiff>searchtol1  && iterM<maxsearchiter  %%%%This loop takes the initial values and
                if abs(tempQMtildemean(2)-tempQ)<=searchtol2  %%%%This is the secondary stopping criterion: if the midpoint T-value is very close to tempT, we take the midpoint, and it doesn't matter how wide the search region is...
                    break %%%%This terminates the while loop
                end
                if tempQMtildemean(2) > tempQ  %%%%If the midpoint gives a T-value above tempT, we can throw away the lower endpoint, which would give a T-value even higher.
                    ThetaL0Mgrid(1)    = ThetaL0Mgrid(2) ;  %%%%Re-define the lower endpoint as the midpoint, and then re-define the midpoint below
                    tempQMtildemean(1) = tempQMtildemean(2) ;

                else  %%%%Otherwise, the midpoint gives a T-value below tempT, and we can throw away the upper endpoint, which would give a T-value even lower.
                    ThetaL0Mgrid(3)    = ThetaL0Mgrid(2) ;  %%%%Re-define the upper endpoint as the midpoint, and then re-define the midpoint below
                    tempQMtildemean(3) = tempQMtildemean(2) ;

                end
                log10Mdiff      = log10(ThetaL0Mgrid(3))-log10(ThetaL0Mgrid(1)) ;
                ThetaL0Mgrid(2) = 10^( log10(ThetaL0Mgrid(1)) + log10Mdiff/2 ) ;
                profitmyopic    = grosspay - ThetaL0Mgrid(2)*commoncost ;
                profnegdiff     = (profitmyopic(2:end,:) - profitmyopic(1:end-1,:)) < 0 ;  %%%%This is a matrix of 1s and 0s indicating whether successive total profits fell from one unit to the next.
                [~,QstopMindx]  = max(profnegdiff) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.
                QstopMindx(QstopMindx==1 & sum(profnegdiff)==0) = length(profnegdiff(:,1));  %%%% EDGE CASE: When a given column of profnedgiff was all zeros,
                                                                                             %%%% meaning profits trended positive throughtout the simulated shock sequence. In this
                                                                                             %%%% case, QstopMindx(ss)==1 and profnegdiff(1,ss)==0.
                tempQMtildemean(2) = mean(QstopMindx-1) ;
                iterM = iterM + 1 ;
            end
            ThetaLlo(ll) = ThetaL0Mgrid(2) ;


            %%%%%%PERFECT FORESIGHT STOP POINTS TO REFINE THE UPPER BOUND OF THE COARSE SEARCH REGION
            for jj=2 %1:length(ThetaL0PFgrid)  %%%%This loop sets the initial values of the bisection search %%%%NOTE: SINCE WE ARE SURE ABOUT THE ENDPOINTS BEING WIDE ENOUGH, WE ONLY NEED THE OBJECTIVE FUNCTION VALUE AT THE MIDPOINT.
                profitmyopic = grosspay - ThetaL0PFgrid(jj)*commoncost ;
                [~,QstopPFindx] = max(profitmyopic) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.
                tempQPFtildemean(jj) = mean(QstopPFindx-1) ;
            end
            while log10PFdiff>searchtol1  && iterPF<maxsearchiter  %%%%This loop takes the initial values and
                if abs(tempQPFtildemean(2)-tempQ)<=searchtol2  %%%%This is the secondary stopping criterion: if the midpoint T-value is very close to tempT, we take the midpoint, and it doesn't matter how wide the search region is...
                    break %%%%This terminates the while loop
                end
                if tempQPFtildemean(2) > tempQ  %%%%If the midpoint gives a T-value above tempT, we can throw away the lower endpoint, which would give a T-value even higher.
                    ThetaL0PFgrid(1)    = ThetaL0PFgrid(2) ;  %%%%Re-define the lower endpoint as the midpoint, and then re-define the midpoint below
                    tempQPFtildemean(1) = tempQPFtildemean(2) ;
                else  %%%%Otherwise, the midpoint gives a T-value below tempT, and we can throw away the upper endpoint, which would give a T-value even lower.
                    ThetaL0PFgrid(3)    = ThetaL0PFgrid(2) ;  %%%%Re-define the upper endpoint as the midpoint, and then re-define the midpoint below
                    tempQPFtildemean(3) = tempQPFtildemean(2) ;
                end
                log10PFdiff      = log10(ThetaL0PFgrid(3))-log10(ThetaL0PFgrid(1)) ;
                ThetaL0PFgrid(2) = 10^( log10(ThetaL0PFgrid(1)) + log10PFdiff/2 ) ;
                profitmyopic     = grosspay - ThetaL0PFgrid(jj)*commoncost ;
                [~,QstopPFindx]  = max(profitmyopic) ; %%%%QstopM now contains the myopic stop point; i.e. the FIRST index after which total profits declined.
                tempQPFtildemean(2) = mean(QstopPFindx-1) ;
                iterPF = iterPF + 1 ;
            end
            ThetaLhi(ll) = ThetaL0PFgrid(2) ;

            tempThetaL = 10.^(logThetaLlowt*log10(ThetaLlo(ll)) + (1-logThetaLlowt)*log10(ThetaLhi(ll))) ;
            ThetaL80(ThetaL80(:,1)==student_id_l(ii),ll+1) = tempThetaL ;
            
            %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
            %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
            %%%%%%%%                           STEP 3, COUNTERFACTUAL SIMULATIONS
            %%%%%%%%  NOTE: for descriptive code comments, see Step 2 above.
            %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
            %%%%%%%%33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333%%%%%%%%
            if tempcon==1
                profitmyopicCF1 = grosspay2(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
                profitmyopicCF2 = grosspay3(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            elseif tempcon==2
                profitmyopicCF1 = grosspay1(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
                profitmyopicCF2 = grosspay3(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            elseif tempcon==3
                profitmyopicCF1 = grosspay1(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
                profitmyopicCF2 = grosspay2(:,Sims3indx) - tempThetaL*commoncost(:,Sims3indx) ;
            end
            
            profnegdiffCF1  = (profitmyopicCF1(2:end,:) - profitmyopicCF1(1:end-1,:)) < 0 ;
            [~,QstopindxCF1] = max(profnegdiffCF1) ;
            QstopindxCF1(QstopindxCF1==1 & sum(profnegdiffCF1)==0) = length(profnegdiffCF1(:,1));
            tempQtildeCF1 = zeros(size(profitmyopicCF1(1,:))) ;
            
            profnegdiffCF2  = (profitmyopicCF2(2:end,:) - profitmyopicCF2(1:end-1,:)) < 0 ;
            [~,QstopindxCF2] = max(profnegdiffCF2) ;
            QstopindxCF2(QstopindxCF2==1 & sum(profnegdiffCF2)==0) = length(profnegdiffCF2(:,1));
            tempQtildeCF2 = zeros(size(profitmyopicCF2(1,:))) ;
            for ss=1:Sims3 
                for kk=1:2  %%%%The value of kk indicates whether we are doing CF1 (the low piece-rate CF) or CF2 (the high piece-rate CF)
                    if kk==1  %%%%This means we are currently doing the low-piece-rate CF
                        Qcurrent = QstopindxCF1(ss)-1 ;
                        if tempcon==1  %%%%In this case, we are doing the low-piece-rate CF which is contract 2, since the actual contract is 1
                            pi0CF = 10 ;
                            pi1CF = 1 ;
                        elseif tempcon>1 %%%%In this case, we are doing the low-piece-rate CF which is contract 1, since the actual contract is either 2 or 3
                            pi0CF = 15 ;
                            pi1CF = 0.75 ;
                        end
                    elseif kk==2  %%%%This means we are currently doing the high-piece-rate CF
                        Qcurrent = QstopindxCF2(ss)-1 ;
                        if tempcon==3  %%%%In this case, we are doing the high-piece-rate CF which is contract 2, since the actual contract is 3
                            pi0CF = 10 ;
                            pi1CF = 1 ;
                        elseif tempcon<3 %%%%In this case, we are doing the high-piece-rate CF which is contract 3, since the actual contract is either 1 or 2
                            pi0CF = 5 ;
                            pi1CF = 1.25 ;
                        end
                    end
                    if Qcurrent<81
                        if Qcurrent==0
                            Tdom3 = Tdom2 ;
                        else
                            Tdom3 = SimTime(Qcurrent,Sims3indx(ss)) + Tdom2 ;
                        end
                        CRev = interp1(Tdom,EQgivenT(:,Qcurrent+1),Tdom2,'linear') ;
                        if Qcurrent==0
                            CRev(CRev<2)  = (CRev(CRev<2)/2)*(2*pi1CF + pi0CF) ;
                            CRev(CRev>=2) = CRev(CRev>=2)*pi1CF + pi0CF ;
                        elseif Qcurrent==1
                            CRev(CRev<1)  = CRev(CRev<1)*(2*pi1CF + pi0CF) ;
                            CRev(CRev>=1) = CRev(CRev>=1)*pi1CF + pi1CF + pi0CF ;
                        else
                            CRev = CRev*pi1CF ;
                        end
                        if max(Tdom3)<=knotvec_c(end)
                            CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qcurrent+1,Sims3indx(ss)) ;
                        else
                            CCost                   = zeros(size(Tdom3)) ;
                            interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                            CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                            CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +...
                                (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                            CCost                   = CCost - commoncost(Qcurrent+1,Sims3indx(ss)) ;
                        end
                        CValcurrent   = max(0,max(CRev - tempThetaL*CCost)) ;
                        if kk==1
                            ValFuncurrent = profitmyopicCF1(Qcurrent+1,ss) + CValcurrent ;
                        elseif kk==2
                            ValFuncurrent = profitmyopicCF2(Qcurrent+1,ss) + CValcurrent ;
                        end
                        stopflag = 0 ;
                        while stopflag==0
                            if CValcurrent==0
                                break
                            end
                            Qnext = Qcurrent + 1 ;
                            Tdom3 = SimTime(Qnext,Sims3indx(ss)) + Tdom2 ;
                            CRev  = interp1(Tdom,EQgivenT(:,Qnext+1),Tdom2,'linear') ; %%%%START WITH EXPECTED OUTPUT q FOR EACH VALUE OF t IN Tdom3
                            if Qnext==1
                                CRev(CRev<1)  = CRev(CRev<1)*(2*pi1CF + pi0CF) ;
                                CRev(CRev>=1) = CRev(CRev>=1)*pi1CF + pi1CF + pi0CF ;
                            else
                                CRev = CRev*pi1CF ;
                            end
                            if max(Tdom3)<=knotvec_c(end)
                                CCost = BsplineEval3(knotvec_c,pi_c,Tdom3) - commoncost(Qnext+1,Sims3indx(ss)) ;
                            else
                                CCost                   = zeros(size(Tdom3)) ;
                                interpindx              = find(Tdom3<=knotvec_c(end),1,'last') ;
                                CCost(1:interpindx)     = BsplineEval3(knotvec_c,pi_c,Tdom3(1:interpindx)) ;
                                CCost(interpindx+1:end) = Cend + Cprimeend*(Tdom3(interpindx+1:end)-knotvec_c(end)) +...
                                    (max(0,Cdblprimeend)/2)*(Tdom3(interpindx+1:end)-knotvec_c(end)).^2 ;
                                CCost                   = CCost  - commoncost(Qcurrent+1,Sims3indx(ss)) ;
                            end
                            CValnext   = max(0,max(CRev - tempThetaL*CCost)) ;
                            if kk==1
                                ValFunnext = profitmyopicCF1(Qnext+1,ss) + CValnext ;
                            elseif kk==2
                                ValFunnext = profitmyopicCF2(Qnext+1,ss) + CValnext ;
                            end
                            if  ValFuncurrent > ValFunnext
                                stopflag = 1 ;
                            else
                                Qcurrent      = Qnext ;
                                CValcurrent   = CValnext ;
                                ValFuncurrent = ValFunnext ;
                            end
                            if Qcurrent+1>=81
                                stopflag = 1 ;
                            end
                        end
                        
                        if kk==1
                            tempQtildeCF1(ss) = Qcurrent ;
                        elseif kk==2
                            tempQtildeCF2(ss) = Qcurrent ;
                        end
                    else
                        if kk==1
                            tempQtildeCF1(ss) = Qcurrent ;
                        elseif kk==2
                            tempQtildeCF2(ss) = Qcurrent ;
                        end
                    end
                end
            end
            
            if tempcon==1
                CFQsamp80_1_2(indx80_1_2+1:indx80_1_2+Sims3) = tempQtildeCF1 ; indx80_1_2 = indx80_1_2 + Sims3 ; 
                CFQsamp80_1_3(indx80_1_3+1:indx80_1_3+Sims3) = tempQtildeCF2 ; indx80_1_3 = indx80_1_3 + Sims3 ; 
            elseif tempcon==2
                CFQsamp80_2_1(indx80_2_1+1:indx80_2_1+Sims3) = tempQtildeCF1 ; indx80_2_1 = indx80_2_1 + Sims3 ; 
                CFQsamp80_2_3(indx80_2_3+1:indx80_2_3+Sims3) = tempQtildeCF2 ; indx80_2_3 = indx80_2_3 + Sims3 ; 
            elseif tempcon==3
                CFQsamp80_3_1(indx80_3_1+1:indx80_3_1+Sims3) = tempQtildeCF1 ; indx80_3_1 = indx80_3_1 + Sims3 ; 
                CFQsamp80_3_2(indx80_3_2+1:indx80_3_2+Sims3) = tempQtildeCF2 ; indx80_3_2 = indx80_3_2 + Sims3 ;
            end
            %%%%IN THE FOLLOWING LOOP, IF NONE OF THE SIMULATED CF STOPPING TIMES ARE BELOW 81, THEN NONE OF THE FOLLOWING
            %%%%ONES WILL BE EITHER, SO WE CAN DIS-CONTINUE THE LOOP FROM ll=1:length(Q80temp)
            if nargout~=2 && min([tempQtildeCF1';tempQtildeCF2'])>=81 && ll<length(Q80temp)
                if tempcon==1
                    CFQsamp80_1_2(indx80_1_2+1:indx80_1_2+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_1_2 = indx80_1_2 + Sims3*(length(Q80temp)-ll) ;
                    CFQsamp80_1_3(indx80_1_3+1:indx80_1_3+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_1_3 = indx80_1_3 + Sims3*(length(Q80temp)-ll) ;
                elseif tempcon==2
                    CFQsamp80_2_1(indx80_2_1+1:indx80_2_1+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_2_1 = indx80_2_1 + Sims3*(length(Q80temp)-ll) ;
                    CFQsamp80_2_3(indx80_2_3+1:indx80_2_3+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_2_3 = indx80_2_3 + Sims3*(length(Q80temp)-ll) ;
                elseif tempcon==3
                    CFQsamp80_3_1(indx80_3_1+1:indx80_3_1+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_3_1 = indx80_3_1 + Sims3*(length(Q80temp)-ll) ;
                    CFQsamp80_3_2(indx80_3_2+1:indx80_3_2+Sims3*(length(Q80temp)-ll)) = 81 ; indx80_3_2 = indx80_3_2 + Sims3*(length(Q80temp)-ll) ;
                end
                break
            end
        end
        
 
    end 

end 

%%
CFQsampACT_1_2(isnan(CFQsampACT_1_2)) = [] ; CFQsampACT_1_2(CFQsampACT_1_2>81) = 81 ; CFQsampACT_1_2(CFQsampACT_1_2<1) = 1 ; 
CFQsampACT_1_3(isnan(CFQsampACT_1_3)) = [] ; CFQsampACT_1_3(CFQsampACT_1_3>81) = 81 ; CFQsampACT_1_3(CFQsampACT_1_3<1) = 1 ; 
CFQsampACT_2_1(isnan(CFQsampACT_2_1)) = [] ; CFQsampACT_2_1(CFQsampACT_2_1>81) = 81 ; CFQsampACT_2_1(CFQsampACT_2_1<1) = 1 ; 
CFQsampACT_2_3(isnan(CFQsampACT_2_3)) = [] ; CFQsampACT_2_3(CFQsampACT_2_3>81) = 81 ; CFQsampACT_2_3(CFQsampACT_2_3<1) = 1 ; 
CFQsampACT_3_1(isnan(CFQsampACT_3_1)) = [] ; CFQsampACT_3_1(CFQsampACT_3_1>81) = 81 ; CFQsampACT_3_1(CFQsampACT_3_1<1) = 1 ; 
CFQsampACT_3_2(isnan(CFQsampACT_3_2)) = [] ; CFQsampACT_3_2(CFQsampACT_3_2>81) = 81 ; CFQsampACT_3_2(CFQsampACT_3_2<1) = 1 ;  

CFQsampMRG_1_2(isnan(CFQsampMRG_1_2)) = [] ; CFQsampMRG_1_2(CFQsampMRG_1_2>81) = 81 ; CFQsampMRG_1_2(CFQsampMRG_1_2<1) = 1 ; 
CFQsampMRG_1_3(isnan(CFQsampMRG_1_3)) = [] ; CFQsampMRG_1_3(CFQsampMRG_1_3>81) = 81 ; CFQsampMRG_1_3(CFQsampMRG_1_3<1) = 1 ; 
CFQsampMRG_2_1(isnan(CFQsampMRG_2_1)) = [] ; CFQsampMRG_2_1(CFQsampMRG_2_1>81) = 81 ; CFQsampMRG_2_1(CFQsampMRG_2_1<1) = 1 ; 
CFQsampMRG_2_3(isnan(CFQsampMRG_2_3)) = [] ; CFQsampMRG_2_3(CFQsampMRG_2_3>81) = 81 ; CFQsampMRG_2_3(CFQsampMRG_2_3<1) = 1 ;   

CFQsamp80_1_2(isnan(CFQsamp80_1_2)) = [] ; CFQsamp80_1_2(CFQsamp80_1_2>81) = 81 ; CFQsamp80_1_2(CFQsamp80_1_2<1) = 1 ; 
CFQsamp80_1_3(isnan(CFQsamp80_1_3)) = [] ; CFQsamp80_1_3(CFQsamp80_1_3>81) = 81 ; CFQsamp80_1_3(CFQsamp80_1_3<1) = 1 ; 
CFQsamp80_2_1(isnan(CFQsamp80_2_1)) = [] ; CFQsamp80_2_1(CFQsamp80_2_1>81) = 81 ; CFQsamp80_2_1(CFQsamp80_2_1<1) = 1 ; 
CFQsamp80_2_3(isnan(CFQsamp80_2_3)) = [] ; CFQsamp80_2_3(CFQsamp80_2_3>81) = 81 ; CFQsamp80_2_3(CFQsamp80_2_3<1) = 1 ; 
CFQsamp80_3_1(isnan(CFQsamp80_3_1)) = [] ; CFQsamp80_3_1(CFQsamp80_3_1>81) = 81 ; CFQsamp80_3_1(CFQsamp80_3_1<1) = 1 ; 
CFQsamp80_3_2(isnan(CFQsamp80_3_2)) = [] ; CFQsamp80_3_2(CFQsamp80_3_2>81) = 81 ; CFQsamp80_3_2(CFQsamp80_3_2<1) = 1 ;  


%%%%THE FOLLOWING ARE WEIGHTS ON EACH OBSERVATION FROM CFQsamp80, TAKING INTO ACCOUNT THAT MULTIPLE TYPES WERE SIMULATED FOR
%%%%EACH PERSON WITH FULL OUTPUT.
weight80_1 = 1/sum(~isnan(Q80freq(:,1,1))) ; %%%%weights for those assigned to contract 1
weight80_2 = 1/sum(~isnan(Q80freq(:,1,2))) ; %%%%weights for those assigned to contract 2
weight80_3 = 1/sum(~isnan(Q80freq(:,1,3))) ; %%%%weights for those assigned to contract 3
 
NActive1 = sum(contract==1 & Q>=2) ; 
NActive2 = sum(contract==2 & Q>=2) ; 
NActive3 = sum(contract==3 & Q>=2) ; 
 
weightMRG_1 = (NActive3-NActive1)/NActive3 ;  %%%%This is the fraction of observations that are "missing" from the active set in group 1
weightMRG_1 = weightMRG_1/(NMarginal1/NActive1) ;  %%%%This now represents how much each head should count for in CFQsampMRG_1_#
weightMRG_2 = (NActive3-NActive2)/NActive3 ;  %%%%This is the fraction of observations that are "missing" from the active set in group 1
weightMRG_2 = weightMRG_2/(NMarginal2/NActive2) ;  %%%%This now represents how much each head should count for in CFQsampMRG_1_#
 
Qdom = (1:81)' ; 
G1mod = zeros(size(Qdom)) ; 
G2mod = zeros(size(Qdom)) ; 
G3mod = zeros(size(Qdom)) ; 
 
for ii=1:length(Qdom)
    G1mod(ii) = weightMRG_2*sum((CFQsampMRG_2_1<=Qdom(ii)))  ...
                  + sum(CFQsampACT_2_1<=Qdom(ii)) + sum(CFQsampACT_3_1<=Qdom(ii)) ...
                  + weight80_2*sum((CFQsamp80_2_1<=Qdom(ii))) + weight80_3*sum((CFQsamp80_3_1<=Qdom(ii))) ;
              
    G2mod(ii) = weightMRG_1*sum((CFQsampMRG_1_2<=Qdom(ii)))  ...
                  + sum(CFQsampACT_1_2<=Qdom(ii)) + sum(CFQsampACT_3_2<=Qdom(ii)) ...
                  + weight80_1*sum((CFQsamp80_1_2<=Qdom(ii))) + weight80_3*sum((CFQsamp80_3_2<=Qdom(ii))) ;
              
    G3mod(ii) = weightMRG_1*sum((CFQsampMRG_1_3<=Qdom(ii))) + weightMRG_2*sum((CFQsampMRG_2_3<=Qdom(ii))) ...
                  + sum(CFQsampACT_1_3<=Qdom(ii)) + sum(CFQsampACT_2_3<=Qdom(ii)) ...
                  + weight80_1*sum((CFQsamp80_1_3<=Qdom(ii))) + weight80_2*sum((CFQsamp80_2_3<=Qdom(ii))) ;
end 
G1mod = G1mod/G1mod(end) ; G1mod = G1mod(2:end-1) ;
G2mod = G2mod/G2mod(end) ; G2mod = G2mod(2:end-1) ;
G3mod = G3mod/G3mod(end) ; G3mod = G3mod(2:end-1) ;
Qdom  = Qdom(2:end-1) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%THIS FINAL PART PUTS ALL OF THE ABOVE INFORMATION TOGETHER
%%%%TO COMPUTE THE VALUE OF THE OBJECTIVE FUNCTION:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Qdom2 = Qdom ;
for ii=1:8
    Qdom2 = sort([Qdom2; (Qdom2(1:end-1) + Qdom2(2:end))/2]) ; 
end

ConfInt = 0.90 ; % 0.95 ; % 0.99 ; % 
CritVal = abs(norminv((1-ConfInt)/2)) ;

Q1samp     = [ones(round(weightMRG_1*sum(Q==1&contract==1)),1); Q(Q>=2 & contract==1)] ; 
[G1emp,Q1] = ecdf(Q1samp) ; Q1 = Q1(2:end); G1emp = G1emp(2:end) ;
if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11
    G1emp80    = max(G1emp(end-1),BsplineEval3(knotvec_q,betaq1,80)) ;
else 
    G1emp80    = interp1([Q1(end-1),max(Q80freq(:,1,1))],[G1emp(end-1),1],80,'linear') ; 
end 
G1emp(end) = G1emp80 ; 
G1lo       = G1emp - CritVal*sqrt( G1emp.*(1-G1emp)/(length(Q1samp)-1) ) ; 
G1up       = G1emp + CritVal*sqrt( G1emp.*(1-G1emp)/(length(Q1samp)-1) ) ;  

Q2samp     = [ones(round(weightMRG_2*sum(Q==1&contract==2)),1); Q(Q>=2 & contract==2)] ; 
[G2emp,Q2] = ecdf(Q2samp) ; Q2 = Q2(2:end); G2emp = G2emp(2:end) ;
if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11 
    G2emp80    = max(G2emp(end-1),BsplineEval3(knotvec_q,betaq2,80)) ;
else 
    G2emp80    = interp1([Q2(end-1),max(Q80freq(:,1,2))],[G2emp(end-1),1],80,'linear') ; 
end 
G2emp(end) = G2emp80 ; 
G2lo       = G2emp - CritVal*sqrt( G2emp.*(1-G2emp)/(length(Q2samp)-1) ) ; 
G2up       = G2emp + CritVal*sqrt( G2emp.*(1-G2emp)/(length(Q2samp)-1) ) ; 

Q3samp     = Q(Q>=2 & contract==3) ; 
[G3emp,Q3] = ecdf(Q3samp) ; Q3 = Q3(2:end); G3emp = G3emp(2:end) ;
if bootstrap(1)==0 || bootstrap(1)==10 || bootstrap(1)==11 
    G3emp80    = max(G3emp(end-1),BsplineEval3(knotvec_q,betaq3,80)) ;
else 
    G3emp80    = interp1([Q3(end-1),max(Q80freq(:,1,3))],[G3emp(end-1),1],80,'linear') ; 
end 
G3emp(end) = G3emp80 ; 
G3lo       = G3emp - CritVal*sqrt( G3emp.*(1-G3emp)/(length(Q3samp)-1) ) ; 
G3up       = G3emp + CritVal*sqrt( G3emp.*(1-G3emp)/(length(Q3samp)-1) ) ; 




%%



lowtailindx = find(G2emp<=1/3,1,'last') ;
upptailindx = find(G2emp>=2/3,1,'first') ;
ssr_1lowtail = sum( ( (interp1(Qdom,G1mod,Qdom2(1:lowtailindx),'linear') - interp1(Q1,G1emp,Qdom2(1:lowtailindx),'linear') ).^2) ) ;
ssr_2lowtail = sum( ( (interp1(Qdom,G2mod,Qdom2(1:lowtailindx),'linear') - interp1(Q2,G2emp,Qdom2(1:lowtailindx),'linear') ).^2) ) ;
ssr_3lowtail = sum( ( (interp1(Qdom,G3mod,Qdom2(1:lowtailindx),'linear') - interp1(Q3,G3emp,Qdom2(1:lowtailindx),'linear') ).^2) ) ;

ssr_1mid = sum( ( (interp1(Qdom,G1mod,Qdom2(lowtailindx+1:upptailindx-1),'linear') - interp1(Q1,G1emp,Qdom2(lowtailindx+1:upptailindx-1),'linear') ).^2) ) ;
ssr_2mid = sum( ( (interp1(Qdom,G2mod,Qdom2(lowtailindx+1:upptailindx-1),'linear') - interp1(Q2,G2emp,Qdom2(lowtailindx+1:upptailindx-1),'linear') ).^2) ) ;
ssr_3mid = sum( ( (interp1(Qdom,G3mod,Qdom2(lowtailindx+1:upptailindx-1),'linear') - interp1(Q3,G3emp,Qdom2(lowtailindx+1:upptailindx-1),'linear') ).^2) ) ;

ssr_1upptail = sum( ( (interp1(Qdom,G1mod,Qdom2(upptailindx:end),'linear') - interp1(Q1,G1emp,Qdom2(upptailindx:end),'linear') ).^2) ) ;
ssr_2upptail = sum( ( (interp1(Qdom,G2mod,Qdom2(upptailindx:end),'linear') - interp1(Q2,G2emp,Qdom2(upptailindx:end),'linear') ).^2) ) ;
ssr_3upptail = sum( ( (interp1(Qdom,G3mod,Qdom2(upptailindx:end),'linear') - interp1(Q3,G3emp,Qdom2(upptailindx:end),'linear') ).^2) ) ; 


test1up      = interp1(Qdom,G1mod,Qdom2,'linear') - interp1(Q1,G1up,Qdom2,'linear') ; %%%%When G1mod rises above the upper 90% conf band (i.e., when test1 goes positive) we impose a quadratic penalty.
guardrail1up = sum( (test1up.*(test1up>0)).^2 ) ;
test2up      = interp1(Qdom,G2mod,Qdom2,'linear') - interp1(Q2,G2up,Qdom2,'linear') ; %%%%When G2mod rises above the upper 90% conf band (i.e., when test1 goes positive) we impose a quadratic penalty.
guardrail2up = sum( (test2up.*(test2up>0)).^2 ) ;
test3up      = interp1(Qdom,G3mod,Qdom2,'linear') - interp1(Q3,G3up,Qdom2,'linear') ; %%%%When G3mod rises above the upper 90% conf band (i.e., when test1 goes positive) we impose a quadratic penalty.
guardrail3up = sum( (test3up.*(test3up>0)).^2 ) ;

test1lo      = interp1(Qdom,G1mod,Qdom2,'linear') - interp1(Q1,G1lo,Qdom2,'linear') ; %%%%When G1mod rises above the upper 90% conf band (i.e., when test1 goes negative) we impose a quadratic penalty.
guardrail1lo = sum( (test1lo.*(test1lo<0)).^2 ) ;
test2lo      = interp1(Qdom,G2mod,Qdom2,'linear') - interp1(Q2,G2lo,Qdom2,'linear') ; %%%%When G2mod rises above the upper 90% conf band (i.e., when test1 goes negative) we impose a quadratic penalty.
guardrail2lo = sum( (test2lo.*(test2lo<0)).^2 ) ; 
test3lo      = interp1(Qdom,G3mod,Qdom2,'linear') - interp1(Q3,G3lo,Qdom2,'linear') ; %%%%When G3mod rises above the upper 90% conf band (i.e., when test1 goes negative) we impose a quadratic penalty.
guardrail3lo = sum( (test3lo.*(test3lo<0)).^2 ) ;

ssr = ssr_1lowtail       + ssr_2lowtail      + ssr_3lowtail ...
     + ssr_1mid          + ssr_2mid          + ssr_3mid  ...
     + ssr_1upptail      + ssr_2upptail      + ssr_3upptail ...
     + 1000*guardrail1up + 1000*guardrail2up + 1000*guardrail3up... 
     + 1000*guardrail1lo + 1000*guardrail2lo + 1000*guardrail3lo ;  

if nargout==2
    out.ssr_1lowtail = ssr_1lowtail ; 
    out.ssr_2lowtail = ssr_2lowtail ; 
    out.ssr_3lowtail = ssr_3lowtail ; 

    out.ssr_1mid = ssr_1mid ; 
    out.ssr_2mid = ssr_2mid ; 
    out.ssr_3mid = ssr_3mid ; 

    out.ssr_1upptail = ssr_1upptail ; 
    out.ssr_2upptail = ssr_2upptail ; 
    out.ssr_3upptail = ssr_3upptail ; 

    out.guardrail1up = guardrail1up ; 
    out.guardrail2up = guardrail2up ; 
    out.guardrail3up = guardrail3up ; 

    out.guardrail1lo = guardrail1lo ; 
    out.guardrail2lo = guardrail2lo ; 
    out.guardrail3lo = guardrail3lo ; 
    out.ssr = ssr ; 
end



if Display==1

    ConfInt = 0.90 ; % 0.99 ; % 0.95 ; % 
    CritVal = abs(norminv((1-ConfInt)/2)) ;

    Q1samp     = [ones(round(weightMRG_1*sum(Q==1&contract==1)),1); Q(Q>=2 & contract==1)] ;
    [G1emp,Q1] = ecdf(Q1samp) ; Q1 = Q1(2:end); G1emp = G1emp(2:end) ;
    G1emp80    = max(G1emp(end-1),BsplineEval3(knotvec_q,betaq1,80)) ;
    G1emp(end) = G1emp80 ;
    G1lo       = G1emp - CritVal*sqrt( G1emp.*(1-G1emp)/(length(Q1samp)-1) ) ;
    G1up       = G1emp + CritVal*sqrt( G1emp.*(1-G1emp)/(length(Q1samp)-1) ) ;

    Q2samp     = [ones(round(weightMRG_2*sum(Q==1&contract==2)),1); Q(Q>=2 & contract==2)] ;
    [G2emp,Q2] = ecdf(Q2samp) ; Q2 = Q2(2:end); G2emp = G2emp(2:end) ;
    G2emp80    = max(G2emp(end-1),BsplineEval3(knotvec_q,betaq2,80)) ;
    G2emp(end) = G2emp80 ;
    G2lo       = G2emp - CritVal*sqrt( G2emp.*(1-G2emp)/(length(Q2samp)-1) ) ;
    G2up       = G2emp + CritVal*sqrt( G2emp.*(1-G2emp)/(length(Q2samp)-1) ) ;

    Q3samp     = Q(Q>=2 & contract==3) ;
    [G3emp,Q3] = ecdf(Q3samp) ; Q3 = Q3(2:end); G3emp = G3emp(2:end) ;
    G3emp80    = max(G3emp(end-1),BsplineEval3(knotvec_q,betaq3,80)) ;
    G3emp(end) = G3emp80 ;
    G3lo       = G3emp - CritVal*sqrt( G3emp.*(1-G3emp)/(length(Q3samp)-1) ) ;
    G3up       = G3emp + CritVal*sqrt( G3emp.*(1-G3emp)/(length(Q3samp)-1) ) ;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%                FIGURE 6:                  %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure;
    subplot(1,3,1)
    hold on; box on; grid on

    plot(Q1(2:end),G1emp(2:end),'r','LineWidth',3) 
    plot(Q1(2:end),G1lo(2:end),':r','LineWidth',1.5) 
    plot(Q1(2:end),G1up(2:end),':r','LineWidth',1.5) 

    plot(Qdom2,interp1(Qdom,G1mod,Qdom2,'linear'),'--r','LineWidth',2) 
    if ConfInt==0.9
        legend('Empirical','Lower 90% Conf. Bound','Upper 90% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.95
        legend('Empirical','Lower 95% Conf. Bound','Upper 95% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.99
        legend('Empirical','Lower 99% Conf. Bound','Upper 99% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    end
    xlabel('Work Volume A (Contract Group 1)','FontWeight','bold','FontSize',18)
    ylabel('CDFs','FontWeight','bold','FontSize',18)

    subplot(1,3,2)
    hold on; box on; grid on

    plot(Q2(2:end),G2emp(2:end),'b','LineWidth',3) 
    plot(Q2(2:end),G2lo(2:end),':b','LineWidth',1.5) 
    plot(Q2(2:end),G2up(2:end),':b','LineWidth',1.5) 

    plot(Qdom2,interp1(Qdom,G2mod,Qdom2,'linear'),'--b','LineWidth',2)
    if ConfInt==0.9
        legend('Empirical','Lower 90% Conf. Bound','Upper 90% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.95
        legend('Empirical','Lower 95% Conf. Bound','Upper 95% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.99
        legend('Empirical','Lower 99% Conf. Bound','Upper 99% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    end
    xlabel('Work Volume A (Contract Group 2)','FontWeight','bold','FontSize',18)
    ylabel('CDFs','FontWeight','bold','FontSize',18)

    subplot(1,3,3)
    hold on; box on; grid on

    plot(Q3,G3emp,'g','LineWidth',3) 
    plot(Q3,G3lo,':g','LineWidth',1.5) 
    plot(Q3,G3up,':g','LineWidth',1.5) 

    plot(Qdom2,interp1(Qdom,G3mod,Qdom2,'linear'),'--g','LineWidth',2) 
    if ConfInt==0.9
        legend('Empirical','Lower 90% Conf. Bound','Upper 90% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.95
        legend('Empirical','Lower 95% Conf. Bound','Upper 95% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    elseif ConfInt==0.99
        legend('Empirical','Lower 99% Conf. Bound','Upper 99% Conf. Bound','Model','Location','SouthEast','FontWeight','bold','FontSize',16)
    end
    xlabel('Work Volume A (Contract Group 3)','FontWeight','bold','FontSize',18)
    ylabel('CDFs','FontWeight','bold','FontSize',18)

end


%%
if nargout==2
    out.ThetaLACT  = ThetaLACT ;
    out.ThetaLMRG1 = ThetaLMRG1 ;
    out.ThetaLMRG2 = ThetaLMRG2 ;
    out.ThetaL80   = ThetaL80 ;
    out.Q80freq    = Q80freq ;
end





end




